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Abstract 

We apply lattice Boltzmann method to study the phase separation of a two-dimensional 
binary fluid mixture in shear flow. The algorithm can simulate systems described by 
the Navier-Stokes and convection-diffusion equations. We propose a new scheme for 
imposing the shear flow which has the advantage of preserving mass and momentum 
conservation on the boundary walls without introducing slip velocities. Our main re- 
sults concern the presence of two typical lenght scales in the phase separation process, 
corresponding to domains with two different thicknesses. Our simulations at low vis- 
cosity confirm previous results only valid in the limit of infinite viscosity. 



1 Introduction 

When a fluid mixture is suddenly quenched from a disordered initial condition to a coexis- 
tence state below the critical temperature, the two fluids segregate with domains growing 
with time (see, e.g., Gunton et al. 1983). The morphology of these domains is strongly 
influenced by an apphed flow. In the case of a shear flow (for a review see Onuki 1997), the 
domains are observed to grow greatly elongated in the flow direction (see, e.g., Hashimoto 
et al. 1995). The effects of the shear have been considerably investigated in the last years. 
Several numerical simulations conflrm the anisotropic growth of the domains (Ohta et al. 
1990, Rothman 1991, Olson & Rothman 1995, Wu et al. 1997, Padilla & Toxvaerd 1997, 
Shou & Chakrabarti 2000). New recent theoretical results by Corberi et al. (1998)-(1999) 
show that, as an effect of the shear, the mixture contains domains of two different thicknesses 
and that the relative abundance of these domains changes periodically with the logarithm 
of time. However in the papers of Corberi et al. (1998)-(1999) phase separation occurred 
only for the effects of simple diffusion. The main purpose of this paper is to see how the 
hydrodynamics affects this picture. 

The presence of domains with different scales is a very peculiar result for phase separating 
mixtures. Indeed, in the general picture for the unsheared case, a single time-dependent 
length scale i?(t), which measures the typical size of domains, characterizes the behaviour 
of the mixture (see Bray 1994). This length grows with the power law R{t) ~ t", where 
the value of the growth exponent a is strictly related to the physical mechanism operating 
in the segregation process. For example a = 1/3 in a pure diffusive regime. Therefore the 
existence of two characteristic scales makes the phase separation process very interesting 
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from a theoretical point of view and also with relevant practical implications. 

Looking at the effects of hydrodynamics on this problem is not a simple computational 
task. Domains grow fast in the direction of the flow making finite-size effects quite soon 
relevant in simulations (for a recent review, see, Yeomans 1999). In this paper we use 
a lattice Boltzmann scheme (Orlandini et al. 1995, Swift et al. 1996) to simulate the 
convection-diffusion and Navier-Stokes equations for a binary mixture. This method has been 
found very convenient with respect to other numerical algorithms for fluid mixtures (Chen 
& Doolen 1998) since it allows to reach very large time scales. In particular, an advantage of 
the lattice Boltzmann scheme of Orlandini et al. (1995) and Swift et al. (1996), is that the 
correct fluid equilibrium can be imposed by choosing an appropriate free energy. Moreover, 
differently from other lattice gas methods, it allows to control in a separate way the different 
mechanisms of transport tuning independently the fluid viscosity and the diffusivity. The 
mentioned lattice Boltmann scheme has been used with success in several problems of phase 
separation of binary fluid mixtures without shear in two (Osborn et al. 1995) and three 
dimensions (Kendon et al. 1999), and also in the case of complex fluids (Gonnella et al. 
1997). 

A preliminary important part of our work has been to include proper boundary conditions 
for a shear flow in lattice Boltzmann schemes. This problem, considered also by Wagner 
& Yeomans (1998) and Gates et al. (1999), has been studied in many papers on lattice 
Boltzmann schemes for a single fluid (for a review, see, e. g., Chen & Doolen 1998). We 
discuss previous procedures and propose an improved algorithm which allows to strictly 
conserve mass and momentum on the boundary walls without introducing slip velocities. 
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The outline of the paper is as follows. In the next section we describe the fluid model and 
the lattice Boltzmann scheme. In section 3 we describe and compare different approaches 
for implementing shear boundary conditions in the numerical scheme. In section 4 we give 
our results for the phase separation and finally we draw our conclusions. 

2 The lattice Boltzmann method 

Our simulations are based on a lattice Boltzmann scheme developed by Orlandini et al. 
(1995) and Swift et al. (1996). In this scheme the choice of a free energy determines all the 
thermodinamic properties of the fluid. We start from the description of the free energy of the 
fluid mixture; then we show how the lattice Boltzmann equations have been implemented. 

Equilibrium description of the mixture. The free-energy functional generally used in phase 
separation studies of binary mixtures (Bray 1994) is 



where ip is the order parameter which describes the normalized difference in the densities of 
the two fluids. The polynomial terms are related to the bulk properties of the fluid. While 
the parameter b is always positive, the sign of a distinguishes between a disordered and a 
segregated mixture. The case with a > gives a polynomial with one minimum in the origin 
corresponding to the a state with ip — everywhere while two minima are present when 



a < corresponding to the two pure phases with (p — ±^J—a/b. In the following we will 
consider a deep quench (well below the critical value a — 0) into the coexistence region with 
a = — 1 and 6=1. The gradient term is related to the interfacial properties. The equilibrium 




(1) 
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equal to -v2fv and an interfacial width proportional to V2/t (Rowlinson & Widom 1982). 
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The total fluid density n is not explicitly included in (|l|) because we are dealing with an 
incompressible mixture with n constant throughout the volume of the fluid and throughout 
its motion. 

The thermodynamic properties of the fluid follow from the free energy (|ID (Reichl 1980). 
The chemical potential difference between the two fluids is given by 

An = — = aip + bip^ - ftVV- (2) 

The pressure is not a scalar but a tensor Pq,^ since interfaces in the fluid can exert non- 
isotropic forces (Yang et al. 1976). The scalar part Po follows directly from 

= + f - /.^ (vV) - f (Vv.)^ (3) 

where /(v?) is the free-energy density. In order to calculate the pressure tensor Pajj, one has 
to ensure that Pa/3 obeys the condition of mechanical equilibrium daPap = (Evans 1979). 
A suitable choice is 

-Pq/3 = Po^aP + nda'-pdfBip. (4) 

Lattice Boltzmann scheme. We use a square lattice in which each site has eight nearest 
neighbours. The lattice has horizontal and vertical links of unit length and diagonal links 
of length v^- We denote by L the number of lattice sites in each direction. The variables 
of the lattice Boltzmann algorithm are two sets of distribution functions fi{r,t) and gi{r,t), 
defined on each lattice site r at time t. Each of them is associated with a velocity vector 
ej. Defined At as the simulation time step, the quantities SjAt are constrained to be lattice 



vectors: |e.t|At = 1 for the horizontal and vertical directions {i = 1,2,3,4) and |ej|At = 
for the diagonal directions {i = 5, 6, 7, 8). In figure 1 we show a plaquette of the lattice with 
the velocity vectors. Two functions /o(r, t) and go{r,t), corresponding to the distribution 
components that do not propagate (gq = 0), are also taken into account. 

The distribution functions are related to the total density n, to the fluid velocity u and 
to the density difference <f through 

i i i 

They evolve during the time step At according to a single relaxation-time Boltzmann equa- 
tion (Bhatnagar et al. 1954, Chen et al. 1992): 

/,(r + e,At,t + At)-/,(r,t) = -l[/,(r, t) - /f (r, t)], (6) 

r 

r?i(r + e,At,t + At)-(7i(r,t) = - — [g^{v,t)-gl\v,t)i (7) 

where r and are independent relaxation parameters, fi'^{r,t) and gl'^{v^t) are local equi- 
librium distribution functions. 

The time evolution occurs in two steps: a collision and a propagation step. Firstly, the 
distribution functions arriving at the same time on the same site change according to the 
equation (we will write it only for one of the two distributions) 

/ar,t) = /.(r,t)-i[/,(r,t)-/r(r,t)]. (8) 

Then the distribution functions are moved along the lattice directions according to the rule 

/,(r + e,At,t + At) = /;=(r,t). (9) 

Equation (||) comes from combining equations (||) and (|^). 
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The system is initialized with n = 1, (f = and u = everywhere. The choice of a 
constant value of n is justified by the incompressibility of the mixture. Our simulations 
confirm that n stays constant throughout the system during all the time evolution. The 
initial value of corresponds to a symmetric 50:50 mixture with the two fluids completely 
mixed. These quantities are locally conserved in any collision process and, therefore, we 
require that the local equilibrium distribution functions fulfil the equations 

E(/r-/.) = o Y.ft' = n 

i i 

Y.^9?-9i) = ^ Y.9T = ^ (10) 

i i 

Y.ur-fi)^^ = ^ E/re. = nu 

i i 

Following Orlandini et al. (1995) and Swift et al. (1996) the higher moments of the local 
equilibrium distribution functions are defined so that we can obtain continuum equations 
pertinent to a binary fluid mixture. We define 

E /reiaCi^ = Pap + nUaUp , (11) 
i 

J29i''eia = , (12) 

i 

E 9r(iia(iifi = rA/i5«/3 + (pUaUj3 . (13) 
i 

where F is a coefficient related to the mobility of the fluid. We are considering a mixture 



with the two fluids having the same mechanical properties. The constraint (|l^) expresses the 
fact that the two fluids have the same velocity. The local equilibrium distribution functions 
can be expressed as an expansion in terms of the velocity u (Orlandini et al. 1995, Swift et 
al. 1996): 



/f^ = Ai + BjUaCia + ClU^ + DjUaUfseiaCip + G I ^ajjeiaei^ 2 = 1, 2, 3, 4 (14) 

/r = ^11 + BnUaCia + Ciiu^ + DiiUaUpeiaCip + Gn^apeiaCip i = 5, 6, 7, 8 
Similarly 

Qi'^ = Hi + KjUaCia + Jju^ + Q lUo^upeicCip i = 1, 2, 3, 4 (15) 
gl'' = Hii + KjjUaCia + Jiiu^ + QiiUaUpeicCif} z = 5, 6, 7, 8 

The relations (|10|)-(|13|) can be used to fix the coefficients of these expansions. The results 
are reported in the Appendix. 

Continuum equations. It has been shown by Orlandini et al. (1995) and Swift et al. 
(1996) that the above described lattice Boltzmann scheme simulates the continuity, the 
incompressible Navier- Stokes and the convection-diffusion equations 

daUa = , (16) 

dtUa + UpdfsUa = --daPo + l'^'^Ua , (17) 

n 

dt^ + d^{^u^) = rev'— . (18) 

dip 

A Chapman-Enskog expansion (Chapman & Cowling 1970, Swift et al. 1996) of equations 
(I) and (0) to 0(At2) allows to show that Q 



^ = ^\^(At), e = At(r,-l). (19) 



We choose t^= {1 + l/v^)/2 in order to minimize the correction terms of order At^. Thus 

we are left with two free variables r and T which control the kinematic viscosity v and the 

"'^ Actually, spurious terms appear in equations (p^)-(p^), whose effects are shown negligible in the papers 
of Orlandini et al. (1995) and Swift et al. (1996). 



macroscopic mobility TO, respectively. All the simulations were run with n = 1 (this choice 
corresponds to an interfacial width of approximately three lattice spacings), FG = 0.2 and 
units in which At = 1. 

3 Shear boundary conditions 

The next step in implementing the algorithm is to apply a shear flow on the system. We have 
modified previous lattice Boltzmann schemes for a single fluid where sliding walls moving in 
a lattice direction were used to enforce the fluid a given velocity. The flow is assumed to be 
directed along the x-axis (horizontal on the lattice). In this direction we considered periodic 
boundary conditions. Then departing distribution functions on outward-pointing links re- 
enter the lattice via corresponding inward-pointing links on the opposite boundary with no 
constraint on the macroscopic velocity or density. The walls are on the top and the bottom 
of the lattice. Here one has to face two problems: the correct value of the velocity has to be 
imposed to the fluid, and the distribution functions pointing outward have to be managed in 
such a way that mass conservation for both components and momentum conservation for the 
bulk are always guaranteed. In the following we discuss two possible ways for introducing 
the sliding walls 0. 

Equilibrium scheme. In the first scheme that we have considered, the boundary walls 

are placed on the links beyond the upper and the lower sites at a distance equal to half 

different approach based on Lees-Edwards boundary conditions (Lees & Edwards (f972)) has been 
used by Wagner & Yeomans (1999). These boundary conditions, widely used in molecular dynamics, identify 
the point at {x,y) with the point at {x + jLAt,y + L), where 7 is the shear rate. Lattice Boltzmann 
implementation of Lees-Edwards boundary conditions requires also to change the usual definition of the 
local equilibrium distribution functions so that we have preferred to use the more straightforward approach 
with the sliding walls. 
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lattice spacing, in the spirit of some works on shear boundary conditions in lattice Boltzmann 
methods for one fluid (Cornubert et al. 1991, Ladd 1994). We assign to the local equilibrium 
distributions at the top (t) and the bottom (b) of the lattice the values corresponding to 
the velocities Wx^t = 7 — ^ — , Wy^t = 0, w^fi = —1 — ^ — , = 0, being L — 1 the distance 
between the top and the bottom of the lattice (He et al. 1997). This method is expected 
to be accurate for relaxation parameters near unity, since when r = 1 the collision process 
(see equation (^) simply replaces the distribution function with the local equilibrium value. 
For describing the propagation step let us refer to the upper row and to the /j's. After 
a collision, for each site of the upper row, there are three distributions fS,{t), /|(t), f^{t) 
pointing outward the system. The propagation is implemented with the following scheme: 

/7(r + e3At,t + At) = r,{v,t) 

/4(r,t + At) = /2^(r,t) (20) 
/8(r + eiAt,t + At) = r^{v,t) 

This choice can be justified considering reflection of particles against the boundary wall. 
This allows to guarantee the conservation of mass and momentum, that on each site 9 
distribution functions are still defined, and that propagation between adjacent sites on the 
boundary occurs simultaneously to the propagation on all the inner sites, since f^{t) and 
f^{t) are propagated over a distance -\/2 in a time step and /^(i) over a distance equal to 
one. This method, however, as expected and also shown in the following, does not work well 
for arbitrary values of the relaxation parameters. Therefore we heve developed a second way 
for introducing the sliding walls. 
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Collisional scheme. This second method is an improvement of a scheme proposed by Zou 

& He (1997) for one- fluid systems. In this case the boundary walls are placed on the upper 

and lower rows of sites (see, e.g., Noble et al. 1995, Inamuro et al. 1995). We start from the 

propagation step, as it is realized in the original version by Zou & He (1997), and consider 

again the sites at the top of the lattice. After the propagation the functions /o(t), fiif), 

hit), f2i't), feit) and fs(t) are known on each site of the upper row. One uses equations 

L-1 



(|) to determine fjit), fi{t), f^{t) and n. Requiring that the wall velocities w^^t ~ ^" 2 
Wy^t = are imposed to the fluid, we can write 

f lit) + hit) + hit) = n ~ [hit) + hit) + hit) + hit) + hit) + hit)] 

hit) -hit) = n ^ ^ - [hit) - hit) + hit) - hit)] (21) 
hit) + hit) + hit) = hit) + hit) + hit) 



Consistency of equations ([21]) gives 

n = hit) + hit) + hit) + 2 [hit) + hit) + hit)] (22) 

The system of equations (pl[) reduces to two equations with three unknown variables. To 
close the system of equations the bounce-back rule (Lavallee et al. 1991, Cornubert et al. 
1991) is adopted for the distribution functions normal to the boundary. This means that 
the value of f^it) is fixed assigning to it the known value of its outward-pointing counterpart 
hit)- Then one can solve the system of equations ( ^Tl) finding the solutions 

/4(r,t) = hir,t) 

hir,t) = /6(r,t)-i[A(r,t)-/3(r,t)] + in7^ (23) 
hir,t) = /5(r,t) + i[/i(r,t)-/3(r,t)]-in7^ 
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With this choice for the inward-pointing distributions the desidered momentum at the bound- 
ary is achieved. At this point the colhsion step is apphed to all sites, including the boundary 
ones. Unfortunately this scheme does not allow to conserve exactly the total mass due to the 
fact that the distribution functions f^{t), f^it) and fgit) resulting from the collision step are 
not accounted for any more, and this avoids the exact mass conservation (Zou & He 1997). 

We have improved the above scheme in order to overcome the mass conservation prob- 
lem. We consider again the propagation step and the problem of calculating the unknown 
quantities frit), f^{t) and ff,{t). After the propagation step on each upper site one has 
the distribution functions f^{t), f2(t), f^it) and f^it) coming from the nearest sites, 

fo(t — At), which does not propagate, and f^it — At), /2(t — At) and feit — At), all referred 
to the previous time step, which were lost in the version of Zou & He (1997). Mass will be 
conserved if the total density n on each site is equal to the quantity n given by the sum 

h{t, t-At) = /o(t - At) + /5(t - At) + /2(t - At) + /6(t - At) 

+fi{t) + hit) + f2{t) + Mt) + fsit). (24) 



We also require that equations (^) are fulfilled. In order to impose the constraint that on 
all the boundary sites n = n, we have to introduce an extra degree of freedom in the system 
of equations. We have choosed /o(t) since it does not propagate. The solutions of the system 



of equations (|2T| ) , n = h and (0) are 

/o(r, t)=n- [/i(r, t) + ^(r, t)] - 2 [Uv, t) + ^r, t) + Mr, t)] (25) 

and again the (^). In this way it is guaranteed that 

n(r, t) = n(r,t, t — At) 
12 



nuy{r,t) — 

By this procedure, once the system has been initiahzed, the apphcation of the propagation 
and coUision steps goes on preserving mass and momentum conservation and implementing 
the correct velocity values on the boundaries, as it has also been verified numerically. 

Comparison between the two schemes. Here we compare the validity of the two schemes. 
We have studied how the steady state velocity profile is reached and we have measured the 
slip velocity. The shear is applied as usual from the beginning of the phase separation. 
Preliminarly we have checked that the boundary conditions do not introduce any artificial 
discontinuity on the walls. Therefore we have considered the behaviour of the densities n 
and if along vertical sections normal to the fiow direction. In figure 2 (a)-(b) the plots of n 
and for the middle section of the system, are showed for the coUisional scheme. It can be 
seen that the total density n stays constant and equal to 1 all over the system (its deviations 
from 1 are less than 1/10^) and no boundary effects can be observed. The profile of 
with the fiuctuations corresponding to the presence of interfaces, does not show anything 
pathological, too. Similar results have been obtained with the equilibrium scheme. In figure 
3 we show at consecutive times the x-componcnt of the fluid velocity along a vertical line 
in the middle of the system calculated with the coUisional scheme. The velocity profiles are 
independent on the particular vertical section considered. The steady velocity profile is the 
planar Couette fiow (Landau & Lifshitz 1959) with the x-component of the velocity having 
a linear dependence on the y coordinate. The evolution of the velocity profile in figure 3 is 
very similar to that observed in simple fluids (Schlichting 1979, p. 92). The same evolution 
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is obtained with the equihbrium scheme when r = 1. The situation is different for larger 
values of r, when, with the first scheme, the time needed to reach the steady velocity profile 
is much longer. Actually, it happens that for large enough values of r the linear profile is not 
reached before finite-size effects become relevant in the simulations. This can be understood 
from equation (|). Indeed, if r > 1, /f(r,t) takes more time to reach the local equilibrium 
value /f^(r, t) which contains the information about the velocity of walls. 

A measure of the time tr required to reach the steady state velocity profile can be done 
in this way. We say that the steady state is reached if 



where T is a tolerance set to 10~^. In figure 4 we show the results obtained with the collisional 
scheme. The time U required to reach the steady state decreases at increasing values of r. 
The analytic solution of the Navier-Stokes equation for a single fluid subject to a Couette 
flow suggests a similar behaviour. For a single fluid it can be shown that the non-stationary 
part of the velocity profile behaves as a series of exponential terms exp(— vr^m^z/t/L^) with 
m an integer (Schlichting 1979). For this reason we have reported our results in a log- log 
scale, even if in this plot we do not observe a linear behaviour. 

Next we consider the behaviour of the slip velocity. It is defined as the difference between 
the wall velocity and the fluid velocity along the wall itself and should be negligible. The 
slip velocity can be measured by the maximum relative error. At the top boundary it can 
be calculated as 



Er \ux{r, t + At)- M^(r, t) I + \uy{r, t + At) - Uy{r, t) \ 
Er \ux{r,t)\ + \uy{r,t)\ 



< T 



(27) 



Em = max 



Wx,t 



I = 



1,2,. ..,L 



(28) 
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where again Ux{i, L) is the fluid velocity at the i-th site along the x direction at the top 
boundary and Wx,t = l{L — l)/2 is the wall velocity The results of this measure for different 
values of r and for both schemes are presented in flgure 5. For the coUisional scheme 
~ 10~^ over the explored range [0.7, 20] of values of r (figure 5 (a)). In the equilibrium 
scheme has a minimum value at r = 1 and then increases with r (figure 5 (b)) arriving to 
values of order 0.3. The conclusion is that the first scheme gives less convincing results and 
we have preferred to use the coUisional scheme to run our simulations on phase separation. 

4 Simulations of phase separation 

We have studied the dynamics of phase separation on systems with size L — 256, using 
the values written in section 2 for the parameters of the model and varying 7 in the range 
[0.001,0.01] and r in [0.7,5], being mainly interested in the hydrodynamic regime at low 
viscosity. 

Here we report the results obtained for r = 0.7, which corresponds to a viscosity u — 
0.067, and 7 = 0.005. Similar results have been obtained for other values of parameters. 
A sequence of configurations at different values of the shear strain is shown in the left 
column of figure 6. The first step of the time evolution is the formation of well defined 
interfaces. Then, for a while, domains grow isotropically as in the case without shear. At 
7t = 0.4, for example, the effects of the shear can be observed only near the moving walls 
due to the fact that at this time the hnear profile of the velocity is not yet reached and the 
velocity is greater close to the walls (see figure 3). The deformation induced by the fiow can 
be observed everywhere in the system starting from values 7^ ~ 1 . Later on the domains 
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become very stretched and progressively aligned with the flow direction as it can be clearly 
seen at 7^ = 7 and 11. It is also clear, by looking at these two configurations, that the 
system shows inhomogeneities and in particular domains with different thicknesses. 

A quantitative analysis of the length scales present in the system can be done by con- 
sidering the behaviour of the structure factor C(k, t). This quantity is the one of ex- 
perimental interest being accessible through scattering techniques. It corresponds to the 
Fourier transform of the two-points correlation function and can be computed numerically 
as < </?(k, k, t) >, where (/^(k, t) is the Fourier transform of the order parameter (p(r, t) 
and < ... > refers to an average over an ensemble of different initial realizations of the sys- 
tem. The wave vector k in the reciprocal lattice is given by k = —{nxi + Uyj), where i and 

Ij 

j are two unit vectors in the kx and ky directions, respectively, and rix and Uy range from 
— L/2 to L/2 — 1. It could be useful to remind the behaviour of C(k, t) in the case without 
shear. Its shape is that of a volcano with a contour plot similar to that shown in figure 6 at 
jt — 0.4. The radius of the volcano defines the inverse of a characteristic length representing 
the average size of domains at a given time. As time goes by, the height of this volcano 
increases while its support shrinks towards the origin. This means that the characteristic 
length scale associated with the size of domains increases and the system is more and more 
ordered on that scale. The fact that the structure factor is circular refiects the isotropy of 
the system. 

The shear affects in a significative way the evolution of the structure factor. At 7^ = 1 
the shape of the volcano is deformed into an elliptical structure as it can be seen in the 
contour plot of figure 6. Moreover also the profile of the edge of the volcano is deformed and 
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two maxima can be observed distributed around the main diagonal of the kx — ky plane. In 
the further evolution the two broad maxima at 7^ = 1 become divided into two separate foils 
and on each foil two different peaks can be observed. This can be clearly seen at — 7. 
These two foils becomes closer and closer to the axis kx = Q and increasingly more aligned 
along the ky axis. This pattern corresponds to interfaces very elongated in the direction of 
the flow as it appears in the corresponding pictures of configurations. 

The presence of two peaks on each foil implies the existence of different length scales in 
the systems. Due to the symmetry C(k, t) — C(— k, t) the two foils are perfectly symmetric 
and it is sufficient to consider only the two peaks of one foil. Therefore for each direction 
there are two relevant length scales. Focussing on what happens in the direction normal to 
the flow we arrive to the conclusion that there are domains with two characteristic thicknesses 
in the system. This better explains the inhomogeneities observed in the configurations at 
7^ = 7 and 11. 

Also the dynamical evolution of the 4-peaked structure factor is very peculiar. At 7^ = 7 
the peaks which prevail are the ones with a larger value of ky. This corresponds to a 
prevalence of thin domains in the system, as it can perhaps be observed by directly looking 
at the corresponding configuration. Later on the more isotropic peak with the smaller value 
of /Cj^, which has grown faster than the other, becomes higher. This can be observed in figure 
6 at 7^ = 11 in correspondence of a larger abundance of thicker domains. 

The physical picture of what happens in the system is the following. By stretching the 
interfaces the shear induces a stress in the system. This stress increases more and more 
producing a prevalence of thin domains. At a certain point the domains are broken by the 
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flow in points where the stress is larger. This occurs in a cooperative way with a sequence 
of ruptures producing domains which are less stretched. Then the thicker domains, not 
yet broken by the flow, become more numerous and the more isotropic peak of C(k,t) 
dominates as at 7^ = 11. It was not possible to follow longer in our simulations the alternate 
predominance of the two peaks due to inevitable appearance of flnite-size effects. Previous 
analytical results of Corberi et al. (1998) show that in the diffusive regime the alternate 
prevalence of the two peaks continues periodically on a logarithmic time scale. We cannot 
give results on this. Here we mention that in experiments with polymer solutions by Migler et 
al. (1996), after a shallow quench under the critical temperature, segregation is observed with 
4-peaked structure factors similar to those of flgure 6. Moreover, in previous experimental 
papers by Laufer et al. (1973) and Mani et al. (1991)-(1992) always on polymer solutions, 
the description of the dynamics of the network of domains is similar to that given above. 
We believe that our simulations, where for the flrst time the presence of two length scales is 
considered, help to understand these experimental facts. 

We conclude this section giving the results of an exphcit evaluation of the size of the 
domains. This measure is usually obtained from the momenta of the structure factor. Since 
our system is anisotropic we define the average size of domains in the fiow direction as 

. /dkC(k,f) 

^'W=''/dk|MC(k.t) (2") 

and analogously for Ry{t), in the shear direction normal to the fiow. They are one half of 
the wave length corresponding to the characteristic wave vector. In the case without shear 
lattice Boltzmann simulations with our parameters give — Ry ^ t^/^ (Osborn et al. 1995). 
a = 2/3 is the value of the growth exponent typical of the inertial regime which occurs at 
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very low viscosity (Furukawa 1985). From our simulations we extract the behaviour shown 
in figure 7. As in all the other simulations with a shear fiow (Yeomans 1999), it has not 
been possible to extract the value of the exponent a due to the relevance of finite-size effects 
0. Using an argument proposed by Corberi et al. (1999), based on a renormalization group 
approach (Bray 1990), one should expect in the present case, if dynamical scaling for the 
structure factor holds, that the growth in the shear direction obeys the same power law as 
in the case without shear, Ry{t) ~ t^^'^, while in the x-direction Rx{t) ~ 'jtRy(t) ~ 'jt^^^. 
This means that the growth exponent in the flow direction is related to ay, the growth 
exponent in the shear direction, by the relation = ay + 1, where the extra contribution 
1 with respect to the unsheared case comes from the convective term in the convection- 
diffusion equation (|I8D. Our runs, performed on an Alpha Workstation XPIOOO with 1 GB 
of RAM, do not allow to check these predictions. The results of flgure 7 can be more easily 
related to the above discussed dynamical behaviour of the structure factor with the alternate 
dominance of the two peaks which corresponds to the oscillatory pattern of Rx and Ry. We 
also observe that these oscillations make more difficult a possible evaluation of the growth 
exponents. 

5 Conclusions 

We studied phase separation in sheared binary fiuid mixtures at low viscosity. We intro- 
duced proper boundary conditions to impose the shear fiow on the system. These boundary 

conditions strictly conserve mass and momentum and do not introduce any appreciable slip 

•^An analytical determination of the exponent a is only available for the diffusive regime (Corberi et al. 
(1998) , Rapapa & Bray (1999)). It is found that a = 4/3 in the flow direction, and a = 1/3 in the other 
directions as in the case without shear. 
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velocity. Our main results show that the network of domains is characterized by the existence 
of two typical length scales for each direction. Indeed, domains with different thickness are 
clearly visible in the simulations. Structure factors with four peaks confirm these observa- 
tions. During the time evolution there is an alternate dominance of two of these peaks over 
the other couple. This corresponds to a larger abundance first of thin and later of thicker 
domains. Simulations on much larger scales are needed to estabhsh if the oscillations of the 
peaks only characterize an initial transient or if they continue indefinitely, as it is suggested 
by the results of Corberi et al. (1998) vahd for the diffusive regime at infinite viscosity. 
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Appendix 



A suitable choice of the coefficients in the expansions (|l^)-([l5|) consistent with the conditions 



Ao = n-20An, Aj = AAn, An = ^^^^^^^ (30) 

Ti 

Bj = ABji, Bji = — (31) 
a = -y, Cj = 4Cn, Cn = -^ (32) 

Tl 

Dj = ADjj, Djj = - (33) 

o 

Gl^af3 = '^G II,al3-, G II ,xx= ~G II ,yy = — ; G j i ^xy = ~G j j ^yx = —— (34) 

iO O 

Ho = ^-20Hn, Hj = 4Hjj, Hjj = ^ (35) 
Kj = 4Kn, Kjj = ^ (36) 
^0 = -^, Ji = Uii, Jii = -^ (37) 

Qi = ^Qii, Qii = I (38) 

o 
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Figure Captions 



Figure 1: One plaquette of the lattice used in our model with the velocity vectors 

Figure 2: Plots of the total fluid density n (a) and of the density difference (p (b) along 
a vertical section at t = 32, using the coUisional scheme. Simulations were run with 
L = 128, 7 = 0.01 and r = 1. 

Figure 3: Plot of the x-componcnt of the fluid velocity as a function of the y coordinate at 
consecutive times: (A) t = 2, (o) t — 8, (★) t — 16, (•) t — 32 for the coUisional sheme 
of shear boundary conditions. 

Figure 4: The time tr required to reach the steady state is shown as a function of the 
relaxation parameter r, using the coUisional scheme. 

Figure 5: The maximum relative error in the velocities at the top boundary is reported 
for the coUisional (a) and the equilibrium (b) scheme as a function of r. Measures were 
taken at the times shown in figure 4. 

Figure 6: Configurations of the system are shown in the left column at different values of 
the shear strain jt. In the right column the structure factor is contour-plotted at the 
same values of jt. 

Figure 7: Evolution of the average domain size in the shear (lower curve) and fiow (upper 
curve) directions. The slope oi Rx is 1.1. 
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